{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "from MalardClient.MalardClient import MalardClient\n",
    "from MalardClient.DataSet import DataSet\n",
    "from MalardClient.BoundingBox import BoundingBox\n",
    "from MalardClient.MaskFilter import MaskFilter\n",
    "\n",
    "\n",
    "from osgeo import ogr, osr\n",
    "from datetime import datetime\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "minX=-300000 \n",
    "maxX=-200000\n",
    "minY=-2400000\n",
    "maxY=-2300000\n",
    "minT=1298394420\n",
    "maxT=1298871335 \n",
    "\n",
    "client = MalardClient()\n",
    "ds = DataSet( \"cryotempo\", \"GRIS_BaselineC_Q2\", \"greenland\"  )\n",
    "\n",
    "proj4 = client.getProjection(ds).proj4\n",
    "\n",
    "shapeFile = \"/data/puma1/scratch/cryotempo/masks/icesheets.shp\"\n",
    "\n",
    "drv = ogr.GetDriverByName('ESRI Shapefile') #We will load a shape file\n",
    "ds_in = drv.Open(shapeFile)    #Get the contents of the shape file\n",
    "\n",
    "lyr_in = ds_in.GetLayer(0)    #Get the shape file's first layer\n",
    "\n",
    "#print(lyr_in.GetFeatureCount())\n",
    "outdriver=ogr.GetDriverByName('MEMORY')\n",
    "source=outdriver.CreateDataSource('memData')\n",
    "\n",
    "#open the memory datasource with write access\n",
    "tmp=outdriver.Open('memData',1)\n",
    "\n",
    "# Create a Polygon from the extent tuple\n",
    "ring = ogr.Geometry(ogr.wkbLinearRing)\n",
    "ring.AddPoint(minX,minY)\n",
    "ring.AddPoint(minX, maxY)\n",
    "ring.AddPoint(maxX, maxY)\n",
    "ring.AddPoint(maxX, minY)\n",
    "ring.AddPoint(minX,minY)\n",
    "poly = ogr.Geometry(ogr.wkbPolygon)\n",
    "poly.AddGeometry(ring)\n",
    "\n",
    "outLayer = source.CreateLayer(\"filter\")\n",
    "\n",
    "# Add an ID field\n",
    "idField = ogr.FieldDefn(\"id\", ogr.OFTInteger)\n",
    "outLayer.CreateField(idField)\n",
    "\n",
    "# Create the feature and set values\n",
    "featureDefn = outLayer.GetLayerDefn()\n",
    "feature = ogr.Feature(featureDefn)\n",
    "feature.SetGeometry(poly)\n",
    "feature.SetField(\"id\", 1)\n",
    "outLayer.CreateFeature(feature)\n",
    "feature = None\n",
    "\n",
    "result_layer = source.CreateLayer(\"clippedmask\",geom_type=ogr.wkbPolygon)\n",
    "\n",
    "lyr_in.Clip( outLayer, result_layer  )\n",
    "    \n",
    "def check(x, y, minX, maxX, minY, maxY):\n",
    "    #Create a point\n",
    "    pt = ogr.Geometry(ogr.wkbPoint)\n",
    "    pt.SetPoint_2D(0, x, y)\n",
    "\n",
    "    #Set up a spatial filter such that the only features we see when we\n",
    "    #loop through \"lyr_in\" are those which overlap the point defined above\n",
    "    \n",
    "    result_layer.SetSpatialFilter(pt)\n",
    "    #Loop through the overlapped features and display the field of interest\n",
    "    if result_layer.GetFeatureCount() > 0:\n",
    "        return True\n",
    "    else:\n",
    "        return False\n",
    "    \n",
    "def applyMask( xSrs, ySrs, minX, maxX, minY, maxY ):\n",
    "    withinMask = np.zeros( len(xSrs) )\n",
    "    \n",
    "    for i, xy in enumerate(zip(xSrs, ySrs)):\n",
    "        x,y = xy\n",
    "        withinMask[i] = 1 if check(x,y, minX, maxX, minY, maxY) else 0\n",
    "        \n",
    "    return withinMask\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Started With 212185 points. Finished with 193723.0\n",
      "Masking took 9.033143 in total, 0.146627 to query, 0.108632 to build df, 8.777884 to mask\n"
     ]
    }
   ],
   "source": [
    "#minX, maxX, minY, maxY = lyr_in.GetExtent()\n",
    "\n",
    "bb = BoundingBox(minX, maxX, minY, maxY, minT, maxT )\n",
    "\n",
    "start =datetime.now()\n",
    "resultInfo = client.executeQuery(ds, bb )\n",
    "postQuery = datetime.now()\n",
    "df = resultInfo.to_df\n",
    "postDf = datetime.now()\n",
    "mask = applyMask(df['x'],df['y'],minX,maxX,minY,maxY)\n",
    "postMask = datetime.now()\n",
    "df['withinMask'] = mask\n",
    "\n",
    "print(\"Started With {} points. Finished with {}\".format(len(df),df['withinMask'].sum()))\n",
    "\n",
    "totalDiff = (postMask - start).total_seconds()\n",
    "totalMask = (postMask - postDf).total_seconds()\n",
    "totalDf = (postDf - postQuery).total_seconds()\n",
    "totalQuery = (postQuery - start).total_seconds()\n",
    "\n",
    "print(\"Masking took {} in total, {} to query, {} to build df, {} to mask\".format(totalDiff, totalQuery, totalDf, totalMask))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "193723\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "['File deleted successfully [/data/puma1/scratch/v2/malard/export/cryotempo_GRIS_BaselineC_Q2_-133019006.nc]']"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f = MaskFilter(p_shapeFile=shapeFile)\n",
    "res = client.executeQuery(ds, bb, maskFilters=f )\n",
    "\n",
    "print(len(res.to_df))\n",
    "\n",
    "client.releaseCacheHandle(res.resultFileName)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
